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We use multifractal detrended fluctuation analysis (MF-DFA) , to See query 1 study sunspot number 
fluctuations. The result of the MF-DFA shows that there are three crossover timescales in the fluc- 
tuation function. We discuss how the existence of the crossover timescales is related to a sinusoidal 
trend. Using Fourier detrended fluctuation analysis, the sinusoidal trend is eliminated. The Hurst 
exponent of the time series without the sinusoidal trend is 0.12 ± 0.01. Also we find that these 
fluctuations have multifractal nature. Comparing the MF-DFA results for the remaining data set 
to those for shuffled and surrogate series, we conclude that its multifractal nature is almost entirely 
due to long range correlations. 

Keyboard: New applications of statistical mechanics 



I. INTRODUCTION 



The important feature of the sun's outer regions is the 
existence of a reasonably strong magnetic field. To the 
lowest order of approximation, the sun's magnetic field is 
dipolar in character and is axisymmetric. The strength of 
the field on a typical point on the solar surface is approx- 
imately a few Gauss. There is, however significant vari- 
ation in this value and there are localized regions (called 
sunspots) in which the filed can be much higher [1]. Be- 
cause of the symmetry of the twisted magnetic lines as 
the origin of sunspots. they are generally seen in pairs or 
in groups of pairs at both sides of the solar equator. As 
the sunspot cycle progresses, spots appear closer to the 
sun's equator giving rise to the so called "butterfly dia- 
gram" in the time latitude distribution [2]. The twisted 
magnetic field above sunspots are sites where solar flares 
are observed. It has been found that chromospheric flares 
show a very close statistical relationship with sunspots 
[1]. The number of sunspots is continuously changing in 
time in a random fashion and constitutes a typically ran- 
dom time scries. Figure 1 shows the monthly measured 
number of sunspots in terms of time. The data belongs 
to a data set collected by the Sunspot Index Data Center 
(SIDC) from 1749 up to now [3]. 

Recently, the statistical properties of sun activity have 
been investigated by some methods in chaos theory [4] 
and multifractal analysis [5,6]. The periodical occur- 
rence of hemispheric sunspot numbers have been ana- 
lyzed with respect to the changes in time using wavelet. 
The north-south asymmetries concerning solar activity 
and rotational behavior has been investigated by using 
the wavelet and auto-correlation function [7]. Cross- 
correlation functions between monthly mean sunspot ar- 
eas and sunspot numbers have been determines in some 
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FIG. 1. Observed spot numbers as a function of time. 

papers [8] . The evidence for the existence of " active lon- 
gitudes" on the sun is given by using the autocorrelation 
bmction of daily sunspot numbers [8,9]. See also [10,11] 
about the relation between sunspot number fluctuation 
and number of flares, their evolution step, i.e. duration, 
rise times, decay times, event asymmetries. 

In this paper we would like to characterize the complex 
behavior of sunspot time series through the computation 
of the signal parameters - scaling exponents - which quan- 
tifies the correlation exponents and multifractality of the 
signal. As shown in Figure 1, the sunspot time scries has 
a sinusoidal trend, with a frequency is equal to the well 
known cycle of sun activity, approximately llyears. Be- 
cause of the nonstationary nature of sunspot time series, 
and due to the finiteness of the available data sample, 



we should apply some methods which are insensitive to 
non-stationarities like trends. 

To eliminate the effect of sinusoidal trend, we apply the 
Fourier Dctrcnded Fluctuation Analysis (F-DFA) [12,13]. 
After elimination of the trend we use the Multifractal Dc- 
trcnded Fluctuation Analysis (MF-DFA) to analysis the 
data set. The MF-DFA methods are the modified version 
of detrended fluctuation analysis (DFA) to detect multi- 
fractal properties of time series. The detrended fluctu- 
ation analysis (DFA) method introduced by Peng et al. 
[14] has became a widely-used technique for the deter- 
mination of (mono-) fractal scaling properties and the 
detection of long-range correlations in noisy, nonstation- 
ary time series [14-18]. It has successfully been applied 
to diverse fields such as DNA sequences [14,19], heart 
rate dynamics [20-22], neuron spiking [23], human gait 
[24] , long-time weather records [25] , cloud structure [26] , 
geology [27], ethnology [28], economical time series [29], 
and solid state physics [30] . 

The paper is organized as follows: In Section II we de- 
scribe the MF-DFA and F-DFA methods in detail and 
show that the scaling exponents determined via the MF- 
DFA method are identical to those obtained by the stan- 
dard multifractal formalism based on partition functions. 
We eliminate the sinusoidal trend via the F-DFA tech- 
nique in Section III and investigate the multifractal na- 
ture of the remaining fluctuation. In Section IV, we ex- 
amine the source of multifractality in sunspot data by 
comparison the MF-DFA results for remaining data set 
to those obtained via the MF-DFA for shuffled and sur- 
rogate series. Section V closes with a discussion of the 
present results. 



II. MULTIFRACTAL DETRENDED 
FLUCTUATION ANALYSIS 

The simplest type of the multifractal analysis is based 
upon the standard partition function multifractal for- 
malism, which has been developed for the multifractal 
characterization of normalized, stationary measurements 
[31-34]. Unfortunately, this standard formalism does not 
give correct results for nonstationary time series that are 
affected by trends or that cannot be normalized. Thus, in 
the early 1990s an improved multifractal formalism has 
been developed, the wavelet transform modulus maxima 
(WTMM) method [35], which is based on the wavelet 
analysis and involves tracing the maxima lines in the 
continuous wavelet transform over all scales. The other 
method, the multifractal detrended fluctuation analysis 
(MF-DFA) , is based on the identification of scaling of the 
gth-order moments depending on the signal length and is 
generalization of the standard DFA using only the second 
moment q = 2. 



The MF-DFA does not require the modulus maxima 
procedure in contrast WTMM method, and hence does 
not require more effort in programming and computing 
than the conventional DFA. On the other hand, often 
experimental data are affected by non-stationarities like 
trends, which have to be well distinguished from the in- 
trinsic fluctuations of the system in order to find the cor- 
rect scaling behavior of the fluctuations. In addition very 
often we do not know the reasons for underlying trends 
in collected data and even worse we do not know the 
scales of the underlying trends, also, usually the avail- 
able record data is small. For the reliable detection of 
correlations, it is essential to distinguish trends from the 
fluctuations intrinsic in the data. Hurst rescaled-range 
analysis [36] and other non-detrcnding methods work well 
if the records are long and do not involve trends. But if 
trends are present in the data, they might give wrong 
results. Detrended fluctuation analysis (DFA) is a well- 
established method for determining the scaling behavior 
of noisy data in the presence of trends without knowing 
their origin and shape [14,21,37-39] 

A. Description of the MF-DFA 

The modifled multifractal DFA (MF-DFA) procedure 
consists of five steps. The first three steps are essentially 
identical to the conventional DFA procedure (see e. g. 
[14-18]). Suppose that Xk is a series of length N, and 
that this series is of compact support, i.e. Xk = for an 
insignificant fraction of the values only. 

• Step 1: Determine the "profile" 

i 

Y{i)=Y,[xk-{x)], i^l,...,N. (1) 
fc=i 

Subtraction of the mean {x) is not compulsory, since it 
would be eliminated by the later detrending in the third 
step. 

• Step 2: Divide the profile Y{i) into Ng = int(A^/s) non- 
overlapping segments of equal lengths s. Since the length 
N of the series is often not a multiple of the considered 
time scale s, a short part at the end of the profile may 
remain. In order not to disregard this part of the series, 
the same procedure is repeated starting from the opposite 
end. Thereby, 2Ns segments are obtained altogether. 

• Step 3: Calculate the local trend for each of the 2Ns 
segments by a least-square fit of the series. Then deter- 
mine the variance 

F\s, {y[{^ - 1)^ + - y.i^)? , (2) 

for each segment v, v ~ 1, • ■ • , and 

F\s, i^) ^ - ^ {Y[N -{u- N,)s + ^] - y,{t)f , (3) 
1=1 



for ly = Ns + 1, ■ ■ ■ , 2Ns. Here, is the fitting poly- 

nomial in segment Linear, quadratic, cubic, or higher 
order polynomials can be used in the fitting procedure 
(conventionally called DFAl, DFA2, DFA3, . . .) [14,22]. 
Since the detrending of the time series is done by the 
subtraction of the polynomial fits from the profile, dif- 
ferent order DFA differ in their capability of eliminating 
trends in the series. In (MF-)DFAm [mth order (MF- 
)DFA] trends of order m in the profile (or, equivalently, 
of order m — 1 in the original series) are eliminated. Thus 
a comparison of the results for different orders of DFA al- 
lows one to estimate the type of the polynomial trend in 
the time series [16,17]. 

• Step 4- Average over all segments to obtain the q-th 
order fluctuation function, defined as: 

where, in general, the index variable q can take any real 
value except zero. For q = 2, the standard DFA proce- 
dure is retrieved. Generally we are interested in how the 
generalized q dependent fluctuation functions Fq(s) de- 
pend on the time scale s for different values of q. Hence, 
we must repeat steps 2, 3 and 4 for several time scales 
s. It is apparent that Fq(s) will increase with increasing 
s. Of course, Fg{s) depends on the DFA order m. By 
construction, Fq(s) is only defined for s > jti + 2. 

• Step 5: Determine the scaling behavior of the fluctua- 
tion functions by analyzing log- log plots of Fq{s) versus s 
for each value of q. If the series Xi are long-range power- 
law correlated; Fq{s) increases, for large values of s, as a 
power-law, 

Fqis) ^ s''<^^\ (5) 

In general, the exponent h{q) may depend on q. For 
stationary time series such as fGn (fractional Gaussian 
noise), Y(i) in Eq. 1, will be a fBm (fractional Brownian 
motion) signal, so, < h{q = 2) < 1.0. The exponent 
h[2) is identical to the well-known Hurst exponent H 
[14,15,31]. Also for a nonstationary signal, such as fBm 
noise, Y{i) in Eq. 1, will be a sum of fBm signal, so 
the corresponding scaling exponent of Fg{s) is identified 
by h{q ~ 2) > 1.0 [14,40] (see the appendix for more 
details). In this case the relation between the exponents 
h{2) and H will he H = h{q = 2) - 1. The exponent 
h{q) is known as generalized Hurst exponent. The auto- 
correlation function can be characterized by a power law 
C(s) = (ukUk+s) ~ s~'^ with exponent 7 = 2 — 2H. 
Its power spectra can be characterized by S{uj) ~ o;^^ 
with frequency lo and /3 = 2H — 1, In the nonstationary 
case, correlation exponent and power spectrum scaling 
are 7 = -2H and l3 = 2H + 1, respectively [14,40]. 

For monofractal time series, h(q) is independent of q, 
since the scaling behavior of the variances F'^{s,v) is 



identical for all segments and the averaging procedure 
in Eq. (4) will just give this identical scaling behavior 
for all values of q. If we consider positive values of g, 
the segments v with large variance F'^(s,v) (i. e. large 
deviations from the corresponding fit) will dominate the 
average Fq{s). Thus, for positive values of q, h{q) de- 
scribes the scaling behavior of the segments with large 
fluctuations. For negative values of g, the segments v 
with small variance F'^{s,v) will dominate the average 
Fq{s). Hence, for negative values of g, h{q) describes the 
scaling behavior of the segments with small fluctuations 
[41]. 

B. Relation to standard multifractal analysis 

For a stationary, normalized series the multifractal 
scaling exponents ft.(g) defined in Eq. (5) are directly re- 
lated to the scaling exponents r(g) defined by the stan- 
dard partition function-based multifractal formalism as 
shown below. Suppose that the series Xk of length N is 
a stationary, normalized sequence. Then the detrending 
procedure in step 3 of the MF-DFA method is not re- 
quired, since no trend has to be eliminated. Thus, the 
DFA can be replaced by the standard Fluctuation Anal- 
ysis (FA), which is identical to the DFA except for a 
simplified definition of the variance for each segment 
1/ ^ 1, . . . , Ng. Step 3 now becomes [see Eq. (2)]: 

FMs,iy)^[Y{i.s)-Y{{u-l)s)r. (6) 

Inserting this simplified definition into Eq. (4) and using 
Eq. (5), we obtain 

— ^|r(^s)-r((^-i),s)|4 ^s'^''\ (7) 

For simplicity we can assume that the length N of the 
series is an integer multiple of the scale s, obtaining Ng = 
N/ s and therefore 

N/s 

5]|y(z..)-y((z.-i).)|'^.'^^(')-i. (8) 

i/=i 

This corresponds to the multifractal formalism used e. g. 
in [32,34]. In fact, a hierarchy of exponents Hq simi- 
lar to our h(q) has been introduced based on Eq. (8) in 
[32]. In order to relate also to the standard textbook 
box counting formalism [31,33], we employ the defini- 
tion of the profile in Eq. (1). It is evident that the term 
Y{i's) — Y{{i' — l)s) in Eq. (8) is identical to the sum 
of the numbers Xk within each segment v of size s. This 
sum is known as the box probability ps(j/) in the standard 
multifractal formalism for normalized series Xk, 

us 

pg{v)= Xk^Y{vs)~Y{{v~l)s). (9) 

k={u~l)s+l 



The scaling exponent r(q) is usually defined via the par- 
tition function Zq{s), 



N/s 



(10) 



where g is a real parameter as in the MF-DFA method, 
discussed above. Using Eq. (9) we see that Eq. (10) is 
identical to Eq. (8), and obtain analytically the relation 
between the two sets of multifractal scaling exponents, 



-{q) = qh{q) - 1. 



(11) 



Thus, we observe that h{q) defined in Eq. (5) for the 
MF-DFA is directly related to the classical multifractal 
scaling exponents t((?). Note that h{q) is different from 
the generalized multifractal dimensions 



g-i 



(12) 



that are used instead of r(q) in some papers. While h{q) 
is independent of q for a monofractal time series, D{q) 
depends on q in this case. Another way to characterize a 
multifractal series is the singularity spectrum f{a), that 
is related to T{q) via a Legendre transform [31,33], 



a = t'((j) and f{a)~qa — T{q). 



(13) 



Here, a is the singularity strength or Holder exponent, 
while f{a) denotes the dimension of the subset of the 
series that is characterized by a. Using Eq. (11), we can 
directly relate a and f{a) to h{q), 

a = h{q) + qh'{q) and /(a) = q[a - h{q)] + 1. (14) 

A Holder exponent denotes monofractality, while in 
the multifractal case, the different parts of the structure 
are characterized by different values of a, leading to the 
existence of the spectrum /(a). 



C. Fourier-Detrended Fluctuation Analysis 

In some cases, there exist one or more crossover (time) 
scales Sx separating regimes with different scaling expo- 
nents [16,17]. In this case investigation of the scaling be- 
havior is more complicate and different scaling exponents 
are required for different parts of the series [18]. There- 
fore one needs a multitude of scaling exponents (multi- 
fractality) for a full description of the scaling behavior. 
A crossover usually can arise from a change in the corre- 
lation properties of the signal at different time or space 
scales, or can often arise from trends in the data. To 
remove the crossover due to a trend such as sinusoidal 
trends, Fourier-Detrende Fluctuation Analysis (F-DFA) 
is applied. The F-DFA is a modified approach for the 



analysis of low frequency trends added to a noise in time 
series [12,13,42,43]. 

In order to investigate how we can remove trends hav- 
ing a low frequency periodic behavior, we transform data 
record to Fourier space, then we truncate the first few 
coefficient of the Fourier expansion and inverse Fourier 
transform the series. After removing the sinusoidal 
trends we can obtain the fluctuation exponent by using 
the direct calculation of the MF-DFA. If truncation num- 
bers are sufficient. The crossover due to a sinusoidal trend 
in the log- log plot of Fq{s) versus s disappears. 

III. ANALYSIS OF SUNSPOT TIME SERIES 

As mentioned in section II, a spurious of correlations 
may be detected if time series is nonstationarity, so di- 
rect calculation of correlation behavior, spectral density 
exponent, fractal dimensions etc., don't give the reliable 
results. It can be checked that the sunspot time series 
is nonstationary. One can verified the non-stationarity 
property experimentally by measuring the stability of 
the average and variance in a moving window for exam- 
ple with scale s. Figure 2 shows the standard deviation 
of signal verses scale s isn't saturate. Let us determine 
that whether the data set has a sinusoidal trend or not. 
According to the MF-DFAl method, Generalized Hurst 
exponents /i(g) in Eq. (5) can be found by analyzing log- 
log plots of Fq{s) versus s for each q. Our investigation 
shows that there are three crossover time scales Sy, in 
the log-log plots of Fq{s) versus s for every g's. These 
three crossovers divide Fq{s) into four regions, as shown 
in Figure 3 ( for instance we took q = 2). The existence 
of these regions is due to the competition between noise 
and sinusoidal trend. For s < six and s > sax , the noise 
has the dominating effect [17]. For six < s < S2x and 
S2x < s < S3x, the sinusoidal trend dominates [17]. The 
value of S2x is approximately equal to 130 month which 
is equal to the well known cycle of sun activity. As men- 
tioned before, for very small scales s < six the effect of 
the sinusoidal trend is not pronounced, indicating that 
in this scale region the signal can be considered as noise 
fiuctuating around a constant which is filtered out by 
the MF-DFAl procedure. In this region the generalized 
DFAl exponent is h{q = 2) = 1.12 ± 0.01, where con- 
firms that the process is a non-stationary process with 
anti-correlation behavior. 

To cancel the sinusoidal trend in MF-DFAl, we apply 
F-DFA method for sunspot data. We truncate some of 
the first coefficient of the Fourier expansion of sunspot se- 
ries. According to Figure 4, for eliminating the crossover 
scales, we need to remove approximately 50 terms of the 
Fourier expansion. Then, by inverse Fourier Transforma- 
tion, the noise without sinusoidal trend is extracted. 

The MF-DFAl results of the remanning new signal are 
shown in Figure 5. The sunspot time scries is a multi- 
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FIG. 2. Behavior of standard deviation of sunspot time 
series as a function of time scale. It shows that this time 
series is not stationary and direct calculation of correlation 
gives strongly wrong result. 



fractal process as indicated by the strong q dependence 
of generalized Hurst exponents and r(g) [44]. The q- de- 
pendence of the classical multifractal scaling exponent 
T{q) has different behaviors for q < and g > 0. For 
positive and negative values of q, the slopes of T{q) are 
1.11 ± 0.01 and 1.44 ± 0.01, respectively. According to 
the relation between the Hurst exponent and h(2), i.e. 
h{q = 2) — 1 = H, we find that the Hurst exponent is 
0.12±0.01. This result is equal to value of Hurst exponent 
in small scale of MF-DFAl of noise with sinusoidal trend. 
The fractal dimension is obtained as Df = 2 — H — 1.88 
[16]. The values of derived quantities from MF-DFAl 
method, are given in Table II and Table HI. 

Usually, in the MF-DFA method, deviation from a 
straight line in the log-log plot of Eq. (5) occurs for small 
scales s. This deviation limits the capability of DFA to 
determine the correct correlation behavior for very short 
scales and in the regime of small s. The modified MF- 
DFA is defined as follows [16]: 



mod 



is) 



([Ff"f(s')]2)i/2 5i/2 

([i^f"f(s)]2)l/2 s'1/2 



(fors' > 1), 



(15) 



where ([-Fl'^Hs)]^)^^^ denotes the usual MF-DFA fluctu- 
ation function [defined in Eq. (4)] averaged over several 
configurations of shuffled data taken from the original 
time series, and s' w iV/40. The values of the Hurst 
exponent obtained by modifled MF-DFAl methods for 
sunspot time series is 0.11 ± 0.01. The relative devia- 
tion of the Hurst exponent which is obtained by modified 




FIG. 3. Crossover behavior of log-log plot F{s) versus s 
for sunspot time series for q = 2.0. There are three crossover 
time scales in plot of F{s), at scales six , S2x and ssx . 




FIG. 4. The MF-DFAl functions Fg{s) for sunspot time 
series versus the time scale s in log-log plot. Original time 
series m = 0, truncation of the first 20 terms m = 20 and 50 
terms m — 50. 
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FIG. 5. The 5 dependence of the generahzed Hurst expo- 
nent h{q), the corresponding r((jr) and singularity spectrum 
/(a) are shown in the upper to lower panel respectively for 
sunspot time series without sinusoidal trend. 

MF-DFAl in comparison to MF-DFAl for original data 
is approximately 8.33%. 

TABLE L The values of the Hurst, multifractal scaling and 
generalized multifractal exponents for q — 2.0, for original, 
surrogate and shuffled of temperature fluctuation series ob- 
tained by MF-DFAl. 



Data 


H 


r 


D 


CMB 


0.94 ±0.01 


0.88 ±0.02 


0.88 ±0.02 


Surrogate 


0.88 ±0.01 


0.76 ±0.02 


0.76 ± 0.02 


ShufHcd 


0.50 ±0.001 


0.002 ± 0.002 


0.002 ± 0.002 



IV. COMPARISON OF THE MULTIFRACTALITY 
FOR ORIGINAL, SHUFFLED AND SURROGATE 
SUNSPOT TIME SERIES 

As discussed in the section III the remanning data set 
after the ehmination of sinusoidal trend has the multi- 
fractal nature. In this section we are interested in to 
determine the source of multifractality. In general, two 



different types of multifractality in time scries can be dis- 
tinguished: (i) Multifractality due to a fatness of proba- 
bility density function (PDF) of the time series. In this 
case the multifractality cannot be removed by shuffling 
the series, (ii) Multifractality due to different correla- 
tions in small and large scale fluctuations. In this case 
the data may have a PDF with finite moments, e. g. 
a Gaussian distribution. Thus the corresponding shuf- 
fled time scries will exhibit mono-fractal scaling, since 
all long-range correlations are destroyed by the shuffling 
procedure. If both kinds of multifractality are present, 
the shuffled series will show weaker multifractality than 
the original series. The easiest way to clarify the type of 
multifractality, is by analyzing the corresponding shuf- 
fled and surrogate time series. The shuffling of time se- 
ries destroys the long range correlation. Therefore if the 
multifractality only belongs to the long range correla- 
tion, we should find ft-shufl*?) = 0.5. The multifractality 
nature due to the fatness of the PDF signals is not af- 
fected by the shuffling procedure. On the other hand, 
to determine the multifractality due to the broadness of 
PDF, the phase of discrete fourier transform (DFT) co- 
efficients of sunspot time series are replaced with a set of 
pscudo independent distributed uniform (— tt, tt) quanti- 
ties in the surrogate method. The correlations in the sur- 
rogate series do not change, but the probability function 
changes to the Gaussian distribution. If multifractality 
in the time series is due to a broad PDF, h{q) obtained 
by the surrogate method will be independent of g. If 
both kinds of multifractality are present in sunspot time 
series, the shuffled and surrogate series will show weaker 
multifractality than the original one. 

To check the nature of multifractality, we compare the 
fluctuation function Fq{s), for the original series ( af- 
ter cancelation of sinusoidal trend) with the result of 
the corresponding shufffed, F^^^^{s) and surrogate series 
F™'^{s). Differences between these two fluctuation func- 
tions with the original one, directly indicate the presence 
of long range correlations or broadness of probability den- 
sity function in the original series. These differences can 
be observed in a plot of the ratio Fq{s) / F^^""^ (s) and 
Fq{s) / F^^'^ {s) versus s [44]. Since the anomalous scal- 
ing due to a broad probability density affects Fg{s) and 
F^^^^{s) in the same way, only multifractality due to cor- 
relations will be observed in Fq(s)/F|^"^(s). The scaling 
behavior of these ratios are 

Fq{s)/F^'''{s) - s'*(«)-''-'(«) = s'^PD-'W. (17) 

If only fatness of the PDF is responsible for the multifrac- 
tality, one should find h{q) ~ /ishuf('?) and ft.cor('z) = 0. 
On the other hand, deviations from /icor('?) = indicates 
the presence of correlations, and q dependence of /icoi (?) 
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FIG. 6. Generalized Hurst exponent, h{q) as a function of 
q for original, surrogate and shuffled data. 

indicates that multifractality is due to the long rage cor- 
relation. If only correlation multifractality is present, one 
finds /ishuf (<z) = 0.5. If both distribution and correlation 
multifractality are present, both, /ishuf (?) and ft.sur('z) will 
depend on q. The q dependence of the exponent h{q) for 
original, surrogate and shuffled time series are shown in 
Figures 6. The q dependence of hcor and /ipDF shows that 
the multifractality nature of sunspot time series is due to 
both broadness of the PDF and long range correlation. 
The absolute value of /icorl?) is greater than hp-£)-p{q), 
so the multifractality due to the fatness is weaker than 
the mulifractality due to the correlation. The deviation 
of hsuriq) and /ishuf('7) from h{q) can be determined by 
using test as follows: 

2 ^ Hit) - hoilt)? /-.ON 

2-.aiq,)■^ + a4q^r' ^'^^ 

the symbol " o " can be replaced by " sur" and " shuf ' , 
to determine the confidence level of /igur and /ishuf to 
generalized Hurst exponents of original series, respec- 

2 

tively. The value of reduced chi-square xto = ^ (-^ 
is the number of degree of freedom) for shuffled and sur- 
rogate time series are 1653.47, 1.10, respectively. On 
the other hand the width of singularity spectrum, /(a), 
i.e. Aa = a{qmin) — cxiqmax) for original, surrogate and 
shuffled time series are approximately, 0.44, 0.75 and 0.22 
respectively. These values also show that the multifrac- 
tality due to correlation is dominant [45] . 

The values of the generalized Hurst exponent h{q = 
2.0), multifractal scaling T{q = 2) and generalized mul- 
tifractal exponents {D{q = 2)) for the original, shuffled 
and surrogate of sunspot time series obtained with MF- 
DFAl method are reported in Table II, The related scal- 
ing exponents are indicated in Table III. The values of 
the Hurst exponent obtained by MF-DFAl and modified 



TABLE II. The values of h{q — 2), multifractal scaling 
and generalized multifractal exponents for q — 2.0 for orig- 
inal, surrogate and shuffled of monthly sunspot time series 
obtained by MF-DFAl. 



Data 


h 


r 


D 


Sunspot 


1.12 ±0.01 


1.24 ±0.02 


1.24 ±0.02 


Surrogate 


1.13 ±0.01 


1.26 ±0.02 


1.26 ±0.02 


ShufHed 


0.51 ±0.01 


0.02 ± 0.02 


0.02 ±0.02 



TABLE III. The values of the Hurst (H), power spectrum 
scaling (/3) and auto-correlation scaling (7) exponents for orig- 
inal, surrogate and shuffled of monthly sunspot time series 
obtained by MF-DFAl. 



Data 


H 
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Sunspot 


0.12 ±0.02 


1.24 ±0.02 


-0.24 ±0.02 


Surrogate 


0.13 ±0.02 


1.26 ±0.02 


-0.26 ±0.02 


Shuffled 


0.51 ±0.01 


0.02 ± 0.02 


0.98 ±0.02 



MF-DFAl methods for original, surrogate and shuffied 
sunspot time series are given in Table IV. 

V. CONCLUSION 

The MF-DFA method allows us to determine the mul- 
tifractal characterization of the nonstationary and sta- 
tionary time series. The concept of MF-DFA of sunspot 
time series can be used to gain deeper insight in to the 
processes occurring in nonstationary dynamical system 
such as sunspots formation. We have shown that the 
MF-DFAl result of the monthly sunspot time series has 
three crossover time scale (sx). These crossover time 
scale are due to the sinusoidal trend. To minimizing the 
effect of this trend, we have applied F-DFA on sunspot 
time series. Applying the MF-DFAl method on trun- 
cated data, demonstrated that the monthly sunspot time 
series is a nonstationary time series with anti-correlation 
behavior. The q dependence of h{q) and T(g), indicated 
that the monthly sunspot time series has multifractal be- 
havior. By comparing the generalized Hurst exponent of 
the original time series with the shuffled and surrogate 
one's, we have found that multifractality due to the cor- 
relation has more contribution than the broadness of the 
probability density function. 
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TABLE IV. The value of the Hurst exponent using 
MF-DFAl and modified MF-DFAl for the original, shuffled 
and surrogate of monthly sunspot time series. 



Method 


Sunspot 


Surrogate 


Shuffled 


MF-DFAl 


0.12 ±0.01 


0.13 ±0.01 


0.51 ±0.01 


Modified 


0.11 ±0.01 


0.12 ±0.01 


0.50 ±0.01 



VI. APPENDIX 

In this appendix we derive the relation between the ex- 
ponent h{2) (DFAl exponent) and Hurst exponent of a 
fBm signal. We show that for such nonstationary signal 
the average sample variance (Eq. 4) for q = 2, is propor- 
tional to s''^*\ where h{q = 2) = H + 1. It is shown that 
the averaged sample variance F^{s) behaves as: 



— y\F^{s,iy)] , 



(19) 



where F'^{s,v) is defined as Eq. 2 and Cu is a function 
of Hurst exponent H . 

To prove the statement we note that the data x{k) is a 
fractional Brownian motion (fBm), the partial sums Y{i) 
(Eq. 1) will be a summed fBm signal. In the DFAl, the 
fitting function will have the expression [yi, = a,y + byi). 
The slope hi, and intercept of a least-squares line on 
Y{i) (from to s) for every windows {v) are given by 

s3/12 

i—1 i—1 i—1 



(20) 



respectively. 

Using the Eqs. 4 and 20, the Eq. 19 can be written as 
follows 




12 



Y,iY{z)Y^Y{i 



(21) 



where we have discard the subscript v for simplicity. Now 
let us calculate the functions A, S, C and D in Eq. 21. 
The increment of a summed fBm and fBm signals, i.e. 



x{i) = Y(i)-Y{i-\) 
u(i) ~ x{i) — x{i — 1) 



(22) 



are a fBm x{i) and fGn u{i) noise, respectively. The 
correlation of Y{i) and x{i) are as follows [15] 



^2 



2H-\ 



(H+l) 



H+l 



(23) 



where cr^ = (u{iY). Also the variance of a summed fBm 



signal is (F(i)^) 



j2(H+i) j;^4]^ YiYiWy using the 



Eqs. 20 and 23, it can be easily shown that the Eq. 21 
can be written as follows 



{[F%s^.)])^^Cns'^"+'\ 



(24) 



where C-h is 



4(7^ 



{2H + ?>){H+\f \{H + l){H + 2)f 
12^2 12^2 



[{H+l){H + 2,)f iH + lf{H + 2){H + 3)' 



(25) 



Therefore the standard DFAl exponent for a non- 
stationary signal is related to its Hurst exponent as 
h{q = 2) = H +1. 
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